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Abstract. The idea of this work is to compare a new positive and entropy stable approximate 
Riemann solver by Francois Bouchut with a state-of the-art algorithm for astrophysical fluid 
dynamics. We implemented the new Riemann solver into an astrophysical PPM-code, the 
f^*^ Prometheus code, and also made a version with a different, more theoretically grounded higher 

order algorithm than PPM. We present shock tube tests, two-dimensional instability tests and 
^ forced turbulence simulations in three dimensions. We find subtle differences between the codes 

JjT in the shock tube tests, and in the statistics of the turbulence simulations. The new Riemann 

solver increases the computational speed without significant loss of accuracy. 

1. Introduction 

In modern astrophysics the interplay between observations and numerical experiments plays a 
central role. Typically hydrodynamical flows with high Reynolds numbers and Mach numbers are 
^ studied, and they are modelled by the Euler equations 

-t-t pt + div(pu) = 

{pu)t + div(puu) +Px = pfi 

{pv)t + div(pwu) +py = pf2 

{pw)t + div(pwu) +Pz= pfs 

(1.1) Et + div{{E + p)u) = pi ■ u. 

Here p is the mass density, u = (u, v, w) is the velocity field, p is the pressure, and E is the 



> 



energy density E = \pvi? + pe with e the specific internal energy. External forces are given in 
' units of acceleration by f — (/i, /2, /s)- The system is closed by the equation of state that relates 

^-H p to p and e. In this work we consider ideal gases where pe = p/(7 — 1) for some 7 > 1, and 

isothermal gases. An isothermal gas has constant temperature T, which implies p = a? p for 
a = with R the gas constant and p, the mean molecular weight. In a real astrophysical flow 
J> additional physical phenomena such as magnetic fields, gravitational forces and electromagnetic 

radiation may be important, but this paper is only concerned with the hydrodynamics. The 
specific physical entropy s is defined by the relation 

c3 1 

(1.2) de + pd-=rds 

P 

with T — T{p, e) > the temperature. The second law of thermodynamics implies that 

(1.3) (p0(s))t + div(pu0(s)) <0 

for any smooth, nonincreasing and convex (f>. In high Mach number flows this condition is needed 
to ensure the dissipativity of shocks, since the viscous forces are ignored in ( |1.1| . 



To numerically solve (1.1), shock-capturing finite volume schemes are widely used. In astro- 
physics it is often done with the PPM algorithm described in [4 , often with an iterative method 
to approximate the exact midpoint value of the Riemann fan. This was implemented in the 
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Higher order algorithm: 


Riemann solver: 


Iterative of Prometheus 


HLLC-Bouchut 


PPM reconstruction and back-tracing 


PPM 


PPM-HLLC 


Piecewise hnear in space, Runge-Kutta in time 


RK-exact 


RK-HLLC 



Figure 1.1. The table summarises the four codes we tested. Along the horizontal 
axis the Riemann solver changes, while vertically the higher order algorithm varies 



Prometheus code in 1989, using the iterative Riemann solver of [5] with a fixed number of itera- 
tions, see [5]. An efficient, parallelised version was then implemented in 2001, see [K] and [IB] . 
Stochastic forcing for turbulence simulations was added later, see |TH], [TH]. Results produced by 
the Prometheus code have been presented in many astrophysical publications, for example in |10j . 
|12j . [13] and [15]. We used this code from 2001 as the basis for this work, and implemented our 
changes into it. 

First, we switched the Riemann solver to an HLLC solver with the signal speeds of Bouchut ([T], 
[2]). As far as we know, this is the first implementation of this advancement into an astrophysics 
code. It is well known that approximate Riemann solvers like HLLC are computationally efficient 
and easy to implement. In addition, the Riemann solver of Bouchut has two good properties: 



a) It automatically ensures that a discrete version of the entropy inequality (1.3 1 holds, b) The 
density and pressure stay positive. Often, finite volume codes need to check for each cell update 
whether the new data are physically reasonable, but with these two a priori estimates, no checks 
are necessary apart from underflow treatment. The two estimates are true in the first order case. 

When using this Riemann solver in a higher order scheme, these two properties are not au- 
tomatically inherited. Hence, we introduced a piecewise linear reconstruction, and replaced the 
characteristic backtracing with Runge-Kutta time integration. This was done in such a way that 
positivity is preserved at one half the CFL-number required in the first order case. Such a second- 
order scheme which also satisfies the entropy inequality is however impractical, see [5], but entropy 
stability for first order schemes has so far seemed to be a good condition in practice. A different 
notion of stability comes from scalar conservation laws, which have solutions with nonincreasing 
total variation. This notion is also important in the design of higher order methods for systems. 
The reconstruction algorithm and the Runge-Kutta integration we use ensure a nonincreasing 
total variation when applied to scalar equations. 



This gave rise to four codes as summarized in Figure 1.1 The codes RK-HLLC and PPM-HLLC 
both run about 20% faster than PPM on the same data with the same resolution. However, the 
difference between the algorithms might be larger, since only the original PPM-code was optimised. 
The RK-exact code is the slowest, but it was of fundamental interest in this project to compare 
its accuracy to the RK-HLLC code's. 

In the remainder of this introduction we will first sum up the main ideas of the underlying 
PPM algorithm in the Prometheus code. Then we will describe the new Bouchut-HLLC solver 
and its theoretical advantages. The second order algorithm is outlined next, and we show how it 
preserves positivity. In sections 2,3 and 4 the one-, two- and three-dimensional test comparisons 
are presented. At the end there is a conclusion. 

1.1. PPM and the Prometheus code. The basis of the Prometheus code is the one-dimensional 
PPM- method of 4 with the iterative Riemann solver from \T . Strang splitting is then used to 
handle multidimensions. The crucial point of PPM is the so-called characteristic back-tracing. 
This technique produces a second order approximation to the states at the cell interfaces at the 
half time step, allowing the use of the midpoint method in time. These approximate states are 
then used as input to the Riemann solver. Although the overall accuracy is second order, the 
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spatial reconstruction is piecewise parabolic, which is reported to give better resolution than 
piecewise linear reconstruction. Furthermore, the accuracy at contact discontinuities is improved 
by an algorithm that detects them, and then steepens the reconstructed density. There is also an 
algorithm that adds artificial diffusion in order to avoid oscillations behind strong, slow-moving 
shocks without smearing out the solution much. The reconstruction is required to be monotocity 
preserving. This means that the order of the scheme may drop locally at extremal points of a 
reconstructed quantity, which means a primitive quantity in the case of PPM as in [3] . This drop 
in order is also a feature of the second order reconstruction we use with the Runge-Kutta time 
integration. 

In order to resolve shocks and shock interactions a Riemann problem is solved with the data from 
the back-tracing operation as input. For the Euler equations there is no general explicit formula 
for the solution of the Riemann problem, so in Prometheus an iterative procedure provides instead 
approximate values of the fluxes at cell interfaces. For efficiency the number of iterations is limited 
to a fixed number. 



1.2. The HLLC solver of Bouchut and stability. The notion of an approximate Riemann 
solver goes back to the Roe solver, see [13], which is based on a local linearization of the fluxes 
at the cell interface. We refer to f5T and [5] for a modern presentation of the basic ideas. The 
basic idea is to replace the exact Riemann fans in Godunov's method with something simpler that 
still gives a numerical flux that is consistent and conservative. In addition entropy consistency 
and positivity of density must be somehow ensured. For this linearized solvers require additional 
treatment, a so-called entropy flx. 

The simplest approximate Riemann solver is the HLL solver, where the Riemann fan is replaced 
by a constant state separated by two discontinuities moving with constant speeds Ci and Cr- A 
sufficient condition for stability is that the exact Riemann solution does not have waves with speed 
outside the interval [C/ , Cr] ■ The main weaknesses of this approach is that material contact waves 
are smeared out, and that the signal velocities C/ and Cr have to be guessed. A solution to the 
first problem was hinted at already in 9,, and was carried out in [5D], see also [H], with the so 
called HLLC-solver. The HLLC-solver consists of three discontinuities traveling with speeds C;, 
u* and Cr, where velocity and pressure are held constant across the middle wave, and u* is this 
intermediate value of the velocity. 

In [2] the HLLC-solver was improved by showing that it results from a relaxation system, which 
established its entropy stability. Furthermore, sharp explicit formulas for the signal speeds that 
ensure positivity and entropy stability could be given. We refer to these as Bouchut speeds, and 
they are given by formula (2.133) and Proposition 2.18 of 0. 



1.3. The new second order algorithm. When going to higher order, requiring the reconstruc- 
tion to be entropy dissipative leads to impractical methods, but there is a way to preserve positivity 
in a rigorous manner. In the rewritten version of Prometheus RK-HLLC and RK-exact we used 
the following reconstruction, based on [2j 
As a limiter we use 



(1.4) minmod(a;, a,.) 



0, aittr < 

sign(aj) min (a|a;|, i(|a;| -|- |ar|), a|ar|) , aia^ > 



with a set to 1.4. This is applied to produce the discrete differential 

(1.5) Dpi = ■^minmod(pi+i - pi,pi - Pi-i), 

and Dui, D{u±)i, D{pe)i similarily. The positivity of the reconstructed density is guaranteed since 
Pi — ^ I Dpi I > miufe (pk)- Conservation of momentum dictates that we take the reconstruction slope 

(1.6) D{pu)i = PiDui + uDpi - —DpiDui, 
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and similarily for D(pu±)i. Energy is conserved by replacing with 
(1.7) = e, - ^ (^1 - ^Dpj^ {Duj + {Du^)j) 

when computing the reconstructed internal energy. The extra terms cancel out the conservation 
errors in kinetic energy caused by the linear reconstruction. Hence positivity means that piCi — 
^\D{pe)i\ > 0, or in other words 

(1-8) y {Du^ + {Du^)l) < Ue.. - \\D{pe),\\ -A^ 



To ensure (1.8 1 in a consistent way we first restrict |-D(uj^)i| to less than or equal to A^, and then 
set Duf less than or equal to A^ — |_D(u„)ip. Note that in practice we multiply Aj with a number 



slightly less than one to ensure that the inequality (1.8) is strict. 

We did not apply any special treatment of material contact waves in this code version, and no 
artificial diffusion was added at shocks. 

The numerical time integration is a second order Runge-Kutta method. That is, one does two 
full time steps, and then averages the resulting cell average with the initial one. This procedure 
preserves positivity, and is total variation diminishing. Multidimensionality is taken care of by 
Strang splitting just as in the PPM-codes. 

2. One space dimension: Shock tube tests of Toro 

A basic setup for testing these methods are one-dimensional Ricmann problems, or shock tube 
tests. In five very useful such test problems are given and subjected to several different 
Riemann solvers. The problems are carefully devised to exhibit phenomena known to be hard to 
reproduce numerically. 

As reference solutions we simulated all tests with 10^ grid cells using the original PPM-code. 
In some cases this was not locally converged due to spurious oscillations etc., and we point out 
these anomalies when they occur. In all the runs the CFL-number was 0.4, and we considered 
X S (0, 1) with a resolution of 100 grid cells. 

2.1. Test 1. The first test is not the most severe, but it contains a transsonic rarefaction, which 
noncntropic schemes have trouble with. The initial data are 



(2.1) {p,u,p) 



(1,0.75,1), x < 0.5 

(0.125,0,0.1), a;>0.5. 



All schemes handle the transsonic rarefaction without any signs of a nonentropic glitch, but there 
are differences in the resolution at the rear end of the rarefaction with the PPM doing the best 
job. However PPM gives large oscillations behind the contact discontinuity compared to the other 
codes. With the RK codes there is little difference between the Riemann solvers. We note the 
undershoot in front of the contact, and the less sharp resolution of the contact compared to PPM. 
These observations also hold true at increasing resolution, as illustrated by the L-'^-errors in mass 
density in Figure pT2| 



2.2. Test 2. Test 2 has two rarefactions going apart creating a low density region. The initial 
data are 



(2.2) {p,u,p) 



(1,-2,0.4), a;<0.5 
(1,2,0.4), a; > 0.5 

The solver should be able to handle this without giving negative density or pressure. In particular. 



linearized solvers have trouble with such cases. In the density plots. Figure 2.3 we note a bump in 
the density at a; = 0.5 with the RK-cxact code. We see similar tendencies for the PPM simulation, 
and in the PPM reference solution there is a deep narrow bump. 
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Figure 2.1. Results for Toro test 1. 



Resolution: 


PPM 


PPM-HLLC 


RK-HLLC 


RK-exact 


/i = 0.01 


3.75 


3.11 


5.07 


5.60 


/i = 0.005 


1.48 


1.61 


2.63 


2.92 


h = 0.0025 


0.88 


0.93 


1.73 


1.85 


/i = 0.00125 


0.42 


0.39 


0.95 


0.98 



Figure 2.2. The table shows the L-'^-error in the computed mass density of Toro's 
test 1 for different codes and resolutions. The errors are given in units of 10^'^. 




Figure 2.3. Results for Toro test 2. The results are symmetric around x = 0.5. 



For the RK-HLLC-code positivity was automatically maintained, and it is interesting that we 
get a better approximation of the density value in the middle compared to PPM-HLLC, and also 
of the velocity. Figure 2.4 The front of the rarefaction is however best resolved by the PPM-codes. 

Notice in Figure 2.4 that both PPM and RK-exact (which have the same Riemann solver) 
has oscillations in the velocity near x — 0.5. The RK-exact code especially had problems with 
this test, and positivity had to be artificially imposed for CFL-numbers larger than around 0.05. 
Theoretically a CFL-number less than 0.25 should ensure positivity with an exact Riemann solver, 
so this has to do with the iterative procedure in the Riemann solver not automatically ensuring 
the positivity property. With the iterative solver as part of PPM however, this seemed not to 
cause serious problems. 
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Figure 2.5. Results for Toro test 3. 

2.3. Test 3. Test 3 is a high Mach number shock tube with initial data 

'(1,0,1000), x<0.5 



(2.3) {p,u,p) 



(1,0,0.01), x>0.5. 



We found little difference between the codes here apart from the expected slightly sharper res- 
olution of the PPM-codes, which was most prominent on the contact wave. All codes produced 



spurious effects behind the rarefaction, as seen in the overshoot in the velocity plots, Figure 2.6 



2.4. Test 4. The solution of test 4 has a near stationary shock, that is, the shock speed is small 
compared to the characteristic speeds, hence if the numerical diffusion applied to this shock is 
high, it will be particularly pronounced. The initial data are 



(2.4) ip,u,p) 



(5.99924, 19.5975, 460.894), x < 0.5 
(5.99242, -6.19633, 46.0950), x > 0.5. 



Also, oscillations behind such shocks is a well known phenomenon, and in the original PPM paper, 
a so called 'flattening' algorithm was introduced which essentially adds diffusion. This algorithm 
was used in all the PPM simulations here, but for this test we also tried to switch it off, resulting 
in oscillations in the density of magnitude around 5 percent of the postshock density both for 
PPM and PPM-HLLC. The RK-codes only show small oscillation here, and no special treatment 
was necessary. In the reference PPM solution there are pronounced oscillations both after the 
near stationary shock as well as between the contact and the right moving shock. Both HLLC- 
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Figure 2.7. Results for Toro test 4. 



codes smear the near stationary shock out with 1-2 grid cells more in front of the shock, as seen 



Figure 2.7 (We show the results from PPM with flattening). The difference of 1-2 grid cells was 
maintained when refining to 200 and 400 grid cells. It is probably caused by the signal speeds 
of HLLC-Bouchut slightly overestimating the shock speeds. Otherwise we only note the lower 
resolution of the contact wave with the RK-codes compared to HLLC. 

2.5. Test 5. Test 5 is like test 3 with a background velocity resulting in a near stationary contact 
discontinuity. The initial data are 



(2.5) {p,u,p) 



(1,-19.59745,1000), a; < 0.8 
(1,-19.59745,0.01), a; > 0.8. 



All codes handle this feature reasonably well, which was expected since both Riemann solvers 
exactly resolve contact waves, see Figure |2T8] Note that in this case the RK-codes have comparable 
resolution of the discontinuities to the PPM-codes, but they overshoot the right intermediate 



density slightly. The oscillations occuring behind the rarefaction in the velocity, see Figure 2.9 
is not reported with the 1st order schemes tested by Toro, so it must have something to do 
with the higher order algorithms. It is especially pronounced with PPM- HLLC, but visible in all 
simulations. 
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Figure 2.9. Results for Toro test 5. 

3. Two SPACE DIMENSIONS: MIXING LAYERS 

We now look at transitions from laminar into unstable flows in two dimensions. By the nature of 
the underlying unstable flows, although we observe differences in the output, we are not really able 
to infer much about the quality of the respective codes. However, we observed that the instabilities 
needed to be highly developed before any differences could be seen between the Riemann solvers. 
In other words, at the onset of instability the Riemann solvers seemed to give the same results. 
We recall from the one-dimensional tests that changing the Riemann solver had only a small effect 
on the numerical smearing. This indicates that the 'numerical viscosity' varies little between 
the different schemes, and the sensitivity of our instabilities to numerical viscosity seems to be 
relatively small at their onset. 

3.1. Kelvin-Helmholtz instability. Two layers of fluid moving with different parallel velocities 
are always unstable in the absence of viscosity and external forces. This is referred to as a Kelvin- 
Helmholtz instability, and seems to be an important source of turbulence in many applications. We 
consider a grid-aligned jump in velocity here, and make a small periodic perturbation. The initial 
data are p — 1, j — 1.4, and we let p vary to allow different Mach numbers. The velocity is in the 
y-direction with v = 0.5 for x < 0.5, and v = —0.5 for x > 0.5, however we moved the velocity 
jump one grid cell to the left to break the symmetry. We perturb v with 27re^'^^cos(27ry)/100 
for X < 0.5 and — 27re^^'^^cos(27ry)/100 for x > 0.5. This means we can compute y e (0,1) 
with periodic boundary conditions, and in the x-direction we consider x G (0.1) with reflecting 
boundary conditions. The CFL-number was 0.8 in all simulations. 
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PPM vs. PPM-HLLC RK-exact vs. RK-HLLC 




Figure 3.1. Growth of transversal kinetic energy component. The codes with 
exact solver are represented with dotted lines, and the HLLC version with solid 
lines. We show data from runs with 100^ 200^ and 4002 grid cell, the curve 
steepness increasing with resolution. 



First we take p — I/7, which means that the relative velocity between the layers equals the 
sound speed. We consider the time history of the average of \pv?, as this quantity is often used 

as a measure of the growth of the instability. Figure 
the different codes. We used three different resolutions. 



shows log ^ pv? as a function of time for 
100^, 200^ and 400^ points, and we see 
here that the instability growth rate increases with resolution. This is plausible, since by linear 
instability theory, the growth rate is inverse proportional to the perturbation wavelength. The 
Riemann solver, however, seems to have no influence at all. We measured the CPU-times for the 
different codes for this simulation with 200^ grid points until time t = 2.0. We got the following 
results: 



Code 


PPM 


PPM-HLLC 


RK-HLLC 


RK-exact 


CPU-time 


1.00 


0.81 


0.79 


1.19 



These are the averages of two runs with each code. The numbers are normalised to the result from 
PPM. 



Figure 3.2 shows the time evolution of the velocity field with PPM-HLLC, illustrated by stream- 
lines at times 0.25,0.50 and 1.0. The streamlines were produced with the intrinsic Matlab routine 
'streamslice'. Note that the density of streamlines plotted does not accurately reflect the numer- 
ical resolution, but were chosen to give a clear representation of the observable topological flow 
features. The plots from the original PPM looks very similar. Figure [373] shows the same but this 
time with the RK-HLLC code. Also with these codes we see no significant differences between the 
two Riemann solvers. 

Differences between the schemes only become apparent at later times. Eventually all vortices 
are swallowed by the domain-centered vortex, and we see some differences in at which time ti this 
happens. We used the streamline plots to find approximately when this change in flow topology 
occurs. For example in Figure [3^ we plotted the flow at time t — 20, since the PPM-simulation 
then still had two distinct vortices, while with PPM-HLLC we could only see one. In Figures [375[ - 
3.7 we do the same with different resolutions and codes. Since the resolution and the underlying 
code also influenced the time ii, we chose different plotting times in Figures [3.4|3.7| In all four 
cases considered it is clear that the Riemann solver induces some difference in ti. 
We also considered p 



1 

1007' 



giving a relative Mach number of 10. Figure 



shows the 



growth of the average transversal kinetic energy component. The effect of the Riemann solver is 
indiscernible. Resolution had less influence here than with relative Mach number 1, so we only 
show the data from the 200^-simulation. At Mach numbers this high, the Kelvin-Helmholtz modes 
are linearly stable, and instead of Kelvin-Helmholtz rolls, kink modes develop, see for example IIT] . 
For the Mach 10 case differences are also very subtle. In Figure [3?9l showing filled density contours 
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\ u 






Figure 3.2. The development of the Kelvin- Helmholtz instability is illustrated 
by the streamlines at times 0.25,0.50 and 1.0 from PPM-HLLC with 200^ 
resolution.* 



at time t — 10, one can see slightly more small scale structure with PPM than PPM-HLLC by 



careful inspection of the plots. The similar plots from the RK-codes in Figure 3.10 clearly show 
more smeared out structures than the PPM-codes. There is no noticeable difference between 
RK-HLLC and RK-exact, which means that any effect of changing the Riemann solver is much 
less prominent than the smearing due the RK-algorithm. The superimposed density contours in 
Figure |3.11| also illustrate the increased numerical diffusion of the RK-algorithm, and that this 
suppresses the effect of changing the Riemann solver. 

3.2. Richtmeyer-Meshkov instability. The Richtmeyer-Meshkov instability occurs when a pla- 
nar shock hits a parallel, slightly perturbed density jump. We used the following setup to simulate 
this. The domain was {x,y) € (0,16) x (0,1) with periodic boundary conditions in y, and Neu- 
mann boundary conditions in x. We set up a shock tube problem in the x-direction at x = 1.6 
with density and pressure as in the first Toro test and constant velocity w = — 1 in the x-direction. 
We considered adiabatic gas 7 = 1.4 so as to coincide with the three-dimensional runs rather 
than the shock tube tests. At the line x = f{y) = 3.2 -I- 0.2 cos(27ry), the density fell by a factor 
of 2. When the shock goes through the initial density jump, the boundary f{y) evolves into a 
mushroom-like structure. In Figure [3. 12| we see a slice in x-direction of the density profile at time 
t = 1.0. The shock has just hit the density jump, and we see a weaker shock going through, and 
a reflected wave moving back towards the contact from the shock tube problem. This last wave 
might cause some minor reflected waves, but otherwise the instability is not influenced by other 
features. The boundary at a; = is transparent to the supersonic rarefaction, and we stop before 



*In the plots x is on the horizontal axis, and y on the vertical axis. 
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Figure 3.3. The same as Figure |3^ with RK-HLLC* 




Figure 3.4. Streamhnes from simulations with 200^ cells at time t = 20. PPM- 
HLLC is shown on the left, and PPM to the right.* 

the shock reflection at x = 16 affects the instability. The CFL-number was 0.8. Again it is hard 
to observe differences, see Figure [3.13| but it seems the original PPM resolves the 'extremities' of 
the high density region a bit sharper, at least they extend more. 

We also illustrate the growth of the instability here by the time history of the transversal 
component of kinetic energy in Figure |3.14[ With RK-HLLC the density jump is more smeared 
out before the shock hits it, which explains the less steep slope. 
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Figure 3.5. Streamlines from simulations with 400^ cells a,t t = 18. PPM-HLLC 
is shown on the left, and PPM to the right.* 




Figure 3.6. Streamlines from simulations with200^ cells at time t = 13. RK- 
HLLC is shown on the left, and RK-exact on the right.* 
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Figure 3.8. Growth of transversal kinetic energy component with relative Mach 
number 10. The codes with exact solver are represented with dotted lines, and 
the HLLC version with solid coloured lines. They look the same. 




Figure 3.9. Density contours from Kelvin-Helmholtz instability with Mach num- 
ber 10 at time t = 10. The resolution was 200^ cells. PPM-HLLC is shown on 
the left, and PPM on the right. Density contours range from 0.4 to 1.4.* 




Figure 3.10. The same as in Figure [3^ but with RK-HLLC on the left and 
RK-exact on the right.* 
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Figure 3.11. Superimposed density contours from Kelvin-Hclmholtz instability 
with Mach number 10 at time i = 10 at resolution 200^ cells. On the left PPM is 
represented by green solid lines, and PPM-HLLC with dotted lines. On the right 
RK-exact with green solid lines, and RK-HLLC with dotted lines. RK right.* 




Figure 3.12. Slice of Richtmeyer-Meshkov instability at t = 1, y = 0.49, Com- 
puted with PPM at resolution Ax = Ay = 0.02. 



It is known that the contact wave steepening of PPM may artificially induce instabilities in cer- 
tain cases. For the PPM-codes we observed small scale structures that we believe to be numerical 
noise when repeating the simulation on finer grids. However, by switching off the steepening, we 
got reasonable results. In Figures [3 . 1 Sp . 1 6] we compare versions of PPM and PPM-HLLC without 
steepening. With the highest resolution the density profiles differ strongly, but there is no way to 
tell which code is better. 



The RK-HLLC-code produces a more smeared out structure, see Figure 3.17 



4. Forced isotropic turbulence 

In many real life flows turbulence is an important feature. Since we do not know how to infer 
from simpler test cases how a numerical method will treat turbulence, we now consider simulations 
of actual three-dimensional turbulence. Because of the three-dimensional nature of turbulence, 
to get useful results one needs powerful computational resources, and we were able to perform 
some parallel simulations on the Hitachi SR8000 at the Leibniz Computing Centre in Munich. 
The simulations were part of a larger study on parameters in supersonic turbulence, see |17j . 
We considered the same type of forced isotropic turbulence experiments described in [18]. The 
resolution here was 256'^ equilateral grid cells, and the boundary conditions periodic as in |18j . 
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Figure 3.13. Density contours from Richtmcycr-Mcshkov instability at time t = 
14. We show data from PPM to the left, and PPM-HLLC to the right. Density 
contours range linearly from 0.1 to 0.22. Resolution Ax — Ay — 0.02. 




0.4 0.6 0.8 1 

t 



Figure 3.14. Logarithm of average -^pv^ vs. time t. We see data from PPM 
as a dotted line, PPM-HLLC solid blue and RK-HLLC dashed magenta. The 
resolution was Ax = Ay = 0.02. 




Figure 3.15. The same as in Figure [3.13| but with the contact wave steepening 
turned off in both codes, and resolution Ax = Ay = 0.01. 
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Figure 3.16. The same as in Figure [3.13| but with the contact wave steepening 
turned off in both codes, and resolution Aa; = Ay = 0.005. 
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Figure 3.17. The same as in Figure 3.13 with RK-HLLC. On the left with 



resolution as in Figure 3.13 and on the right with resolution as in Figure 3.15 



We refer to [16] and [18] for details of the experiments and the analysis tools. In addition to \T8\ . 
compressible turbulence simulations with PPM have been investigated by Sytine et al. in [19]. 

The tests consisted of a constant, zero velocity initial state continuously subjected to a stochas- 
tically varying force field f . The forcing was given by evolving its Fourier transform by a so called 
Ornstein-Uhlenbeck process, which is a statistically stationary stochastic process, with parameters 
such that the resulting force was statistically isotropic. Only the larger wavelengths were given a 
nonzero contribution. Note that the Fourier transform of a periodic function can be understood 
as a generalized function given by the coefficients in its Fourier series. By varying the magnitude 
of the forcing, the characteristic velocity of the flow was varied correspondingly. The forcing also 
had a free parameter corresponding to a projection operator regulating the solenoidality of the 
force fleld. For C = 1 the force field is divergence free, and for lower values we have progressively 
stronger compressive force components. We will not study the infiuence of this parameter here, 
just note that all fiows considered were highly compressible. How a gas responds to this injection 
of energy depends a lot on the equation of state, as we will show. 

Since these fiows are highly sensitive to perturbations, it makes no sense to compare the actual 
solutions. Instead we will compare statistical properties of the simulated fiows, since the statistical 
approach has been relatively successful in quantitatively describing turbulence, see for example [7]. 
Note also that each simulation represented a different realisation of the stochastic forcing process. 
One way to extract statistical information is to make a histogram of the different values assumed 
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by a scalar quantity at a fixed time. We can call this to make a probability distribution function 
(PDF). We will consider PDFs for p and the absolute value of the vorticity lo. 

As an indicator of numerical dissipation we will look at the energy spectra, that is, we will look 
at the energy content in each Fourier mode of the velocity field. Parseval's theorem says that the 
total specific kinetic energy equals the integral over the square of the Fourier transformed velocity 
field u(k,t), 

(4.1) / |u(x,t)|2dx = ^u(k,t) • u(k,i)*. 

k 

where •* denotes complex conjugation. In other words, it is given by integrating over the energy 
spectrum function E{k,t), which is defined as the sum of the squares of the Fourier coefficients 
corresponding to each mode where the three-dimensional wave number vector k has absolute value 

k, 

(4.2) E{k,t)^ ^fl(k,i)-ii(k,t)* 

|k|=fc 

times a scaling factor. We refer to [TB] and [TH] for how this was done numerically. 

It is intuitively clear that if the solution has a lot of small scale structure, it indicates low 
numerical diffusion, although spurious oscillations could also play a role. The energy spectrum 
function gives a way to quantify this idea for these highly complex flows, but it is also connected 
to deeper ideas about turbulence, in particular Kolmogorov's theory, see for example [7]. 

Typically a plot of fc i-^ E{k, t) will show three different ranges. For the lowest wave numbers 
the stochastic injection of mechanical energy dominates. Then comes what is called the inertial 
range, where Kolmogorov's theory predicts that E{k,t) drops off as A:~3, due to the transfer of 
energy from vortices of higher to lower length scales. This 'Kolmogorov cascade' has been observed 
for low enough Mach numbers both in experiments and numerical simulations. For the highest 
wave numbers numerical dissipation becomes dominant, and E{k, t) drops off steeply. Between the 
inertial range and the dissipation range, one tends to observe a flattening of E(k,t) in numerical 
simulations. This is called the bottleneck effect and it is still debated whether it has physical 
significance, or whether it is a purely numerical effect, see [B], |18j and the references in these. 
With the resolution here of 256'^ cells, the injection range goes straight into a bottleneck range. 
Since the Kolmogorov theory is derived for incompressible flow, we also deflne the transversal 
energy spectrum Etr{k, t) which only consist of the part of u(k, t) orthogonal to k, so that we only 
take into account the divergence free part of the velocity. 

Some dimensional quantities need to be deflned flrst, but we choose not to go into detail about 
the physical scales as they are not relevant to the code comparisons. Asymptotically the RMS 
(root mean squared) amplitude of the force / would approach {1 — 2( + 3('^)Fo for some prescribed 
value Fq. We use this to deflne the characteristic velocity V by V — {FqL)^^'^, where L is half 
the length I of the sides of the periodic box. Hence V is close to the RMS velocity in the fully 
developed flow. With characteristic Mach number 'Ma' we refer to the ratio of V to the initial 
sound speed. The simulations were run for flve integral time scales T = y. The forcing is strongest 
at wave numbers k such that |k| = fco = and zero for |k| > 2fco- With a we refer to the integer 
1^ = 2, and the initial density is denoted by po- As scaling factor for E{k,t) we take 

The CFL-number was 0.8 in all simulations. 

4.1. Adiabatic gas, characteristic Mach number 17.9, C = 0.1. We flrst compare PPM and 
RK-HLLC on a set-up with an adiabatic equation of state, that is an ideal gas with 7 = 1.4. Most 
of our statistics reveal no signifant difference between the codes, but we see some clear trends in 
the evolution of the energy spectra. The spectra imply that RK-HLLC is more dissipative than 
PPM when the average Mach number is less than about 5. Furthermore the dissipative effects of 
RK-HLLC appears to grow as the Mach number decreases, while the dissipative effects of PPM 
is unaffected by Mach number. 

In the case of an adiabatic gas, the Mach number initially grows sharply, and then falls off 
because the injected kinetic energy dissipates into heat, hence increasing the sound speed, see 
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Figure 4.1. Time history of RMS (root mean squared) Mach number (left) and 
momentum (right) for adiabatic runs. The curve from PPM is labelled 'ad', and 
the curve from RK-HLLC is labelled 'ad (RK)' . 




Figure 4.2. Energy spectrum at the final time t — 5T for the adiabatic runs 
(left), on the right the compensated transversal spectrum. The curves are labelled 
as in Figure |4.1[ 



Figure 4.1 The velocity field behaves statistically as stationary isotropic turbulence after around 



one integral time scale according to Figure 4.3 although even at the termination point t ~ 5T, an 
equilibrium between the energy injection and dissipation was not reached, as that would imply a 
constant RMS momentum in Figure |4.1| 

Figure [42| shows energy spectra at the final time. The energy spectrum function for RK-HLLC 
drops off significantly more sharply for the high wave numbers, and this is to be expected due to 
the less sharp resolution of the RK-HLLC code. Also note the clear bottleneck effect, which is best 
seen in the plot of the ' com pensated' transversal energy spectrum function ^'(fc,t) proportional 



to Etr{k,t)k3 in Figure 4.2 The Kolmogorov theory predicts that should be constant in the 
inertial range, and then drop off in the dissipation range. 

If we look at other times, however, things become more complicated: In Figure [4!3] we see the 
evolution of the energy spectra. 

The fact that in Figure |4.1| the curves differ up to time t = 2T, we attribute to the different 
realisations of the stochastic forcing. Hence one should compare statistics from the different codes 
only for times t > 2T. The RK-HLLC-code gave decreasing energy spectra in time (after t = T), 
while the energy spectra from PPM have no particular time dependency (after t = T). We 
interpret this as an increase in dissipativity of RK-HLLC in time, and we associate it with the 
decrease in Mach number. The observed independence of Mach number for PPM is related to the 
known fact that this scheme retains, and even improves, its accuracy as the advective Courant 
number decreases. This is caused by the high order upwind advection used in PPM. The typical 
advective Courant number in a cell will decrease with the RMS Mach number in these turbulence 
tests, hence PPM should perform well at the lower Mach numbers. 
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Figure 4.3. Time history of transversal energy spectra. Data from PPM are on 
top, and from RK-HLLC underneath. Times are given in units of integral time 
scale T. 



Figure 4.4 shows mass density PDF at different times. There are clear differences between the 
simulations in both the high and low density regions, but this seems to be due to fluctuations 
inherent in the stochastic process behind the forcing, as there is no clear trend. The vorticity 
PDF's in Figure [475| show the same tendency as the energy spectra. From time t = 2T, the tails 
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Figure 4.4. Time history of mass density PDF for adiabatic runs. The curves 
are labelled as in Figure |4.1[ Times are given in units of integral time scale T. 
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Figure 4.5. Time history of \lo\ PDF for adiabatic runs. The curves are labelled 
as in Figure |4.1[ 



in the vorticity PDF's from PPM are clearly longer than those of RK-HLLC, meaning that 
flow contains more small scale vortices. 
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Figure 4.6. Time history of RMS (root mean squared) Mach number and mo- 
mentum for isothermal run with characteristic Mach number 2.1. The curve from 
PPM is labelled 'it\ from PPM-HLLC 'it (HLLC)\ and from RK-HLLC Ht (RK)\ 




Figure 4.7. Time history of RMS (root mean squared) Mach number and mo- 
mentum for isothermal run with characteristic Mach number Ma 21.1. The curves 
are labelled as in Figure [46] 



4.2. Isothermal gas. We also performed simulations with isothermal gas. Here the Mach number 
stays near constant after the initial growth phase, as seen in Figures [OpJl For this reason, we 
only analyse the data from the final time. We found no significant difference between RK-HLLC 
and PPM, and there was even less difference between the two PPM codes. We show PDFs and 
energy spectra from simulations with Mach numbers 2.1 and 21.1 in Figures |4.8|4.10| Again 
bottleneck effects are seen in all simulations, see Figure |4.11[ 

It is a little surprising that the RK-HLLC code does not appear comparably more dissipative 
in these isothermal runs as it did in the adiabatic case, and in the one- and two-dimensional tests. 
For the isothermal test with characteristic Mach number 21.1, it is not so unexpected. This is 
because the RMS Mach number was much higher than in the adiabatic runs, and from Figure 
|4.3| it seems that the difference is less for higher Mach numbers. For the test with characteristic 
Mach number 2.1 however, the RMS Mach number is comparable to that of the adiabatic tests 
at around t — ST. Hence, in addition to Mach number, the equation of state must be taken into 
account when comparing RK-HLLC and PPM. In section [2] we noted that the most significant 
difference between the RK- and PPM-algorithms was that both RK-codes smeared out contact 
discontinuities more. In the isothermal case contact waves are not present, which might explain 
why the codes differ less here than in the adiabatic runs. 



5. Summary 

From our numerical experiments we have made the following observations: 
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Figure 4.8. Density PDFs for isothermal runs from the final time t — ST. On 
the left Ma=2.1, and on the right Ma=21.1. The curves are labelled as in Figure 
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Figure 4.9. Vorticity PDFs for isothermal runs from the final time t — 5T. On 
the left Ma=2.1, and on the right Ma=21.1. The curves are labelled as in Figure 
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Figure 4.10. Energy spectra for isothermal runs from the final time t — 5T. On 
the left Ma=2.1, and on the right Ma=21.1. The curves are labelled as in Figure 

SSI 



• There is slightly more smearing of stationary shocks with HLLC-Bouchut compared to the 
exact solver. 

• The RK-codes smear out most features more than the PPM-codes, and especially contact 
discontinuities. 
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Figure 4.11. Compensated energy spectra for isothermal runs from the final 
time t — 5T. On the left Ma=2.1, and on the right Ma=21.1. The vertical lines 
represent the 'sonic wave number'. Integrating the curve to the left of his lin gives 



cj. The curves are labelled as in Figure 4.6 



• RK-HLLC handles small densities better than the other codes. 

• All codes exhibit spurious oscillations. We see more of them with the PPM-codes, except 
at the near stationary shock, where PPM has a specialised 'flattening' procedure. 

• The growth of Kelvin-Helmholtz and Richtmeyer-Meshkov instabilities appears to be little 
affected by which of the two Riemann solvers are used. 

• In turbulence simulations of adiabatic gas, the dissipativity of RK-HLLC seems to be less 
for higher than for lower Mach numbers, while the dissipation with PPM is independent 
of the Mach number. 

• For turbulence in an adiabatic gas with an RMS (root mean squared) Mach number less 
than about 5, RK-HLLC seems to be more dissipative than PPM. 

• For turbulence with an RMS (root mean squared) Mach number of 2.5 and higher in an 
isothermal gas, there were no significant differences between Riemann solvers or higher 
algorithms. 

The widespread use of PPM in the astrophysics community has lead to concern about how much 
the results depend on this algorithm. We conclude that with respect to the Riemann solver their 
results are accurate. However the efficiency of the HLLC Riemann solver of Bouchut suggests that 
it may be used instead. 
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